Predição da profundidade do solum
Horst
1 de novembro de 2017
rmarkdown::render('main.Rmd', encoding = 'UTF-8', output_dir = "../docs")rm(list = ls())
# Carregar os Pacotes
library(magrittr)## Warning: package 'magrittr' was built under R version 3.4.2library(dplyr)## Warning: package 'dplyr' was built under R version 3.4.2##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(glue)## Warning: package 'glue' was built under R version 3.4.3##
## Attaching package: 'glue'## The following object is masked from 'package:dplyr':
##
## collapselibrary(lattice)
library(latticeExtra)## Warning: package 'latticeExtra' was built under R version 3.4.1## Loading required package: RColorBrewerlibrary(georob)## Warning: package 'georob' was built under R version 3.4.2## Loading required package: parallel## Loading required package: sp## Warning: package 'sp' was built under R version 3.4.1##
## Attaching package: 'georob'## The following object is masked from 'package:stats':
##
## steplibrary(sp)
library(mapview)## Warning: package 'mapview' was built under R version 3.4.1## Loading required package: leaflet## Warning: package 'leaflet' was built under R version 3.4.1library(raster)## Warning: package 'raster' was built under R version 3.4.1##
## Attaching package: 'raster'## The following object is masked from 'package:georob':
##
## cv## The following object is masked from 'package:glue':
##
## trim## The following object is masked from 'package:dplyr':
##
## select## The following object is masked from 'package:magrittr':
##
## extractlibrary(rmarkdown)
library(caret)## Warning: package 'caret' was built under R version 3.4.1## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.4.1##
## Attaching package: 'ggplot2'## The following object is masked from 'package:latticeExtra':
##
## layer# Sistemas de referência de coordenadas (Fonte: http://spatialreference.org/ref/epsg/)
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sirgas2000 <- sp::CRS('+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')
# Rampas de cores
col_soil_var <- topo.colors(100)Caracterização dos dados
Os dados que utilizei para este trabalho são provenientes de um povoamento florestal de 108 ha de \(Pinus taeda\) L. de 29 anos. A área pertence ao município de Campo Belo do Sul, região serrana do Estado de Santa Catarina, Brasil. O clima é do tipo Cfb, mesotérmico, subtropical úmido e com precipitação média de 1.647 mm com chuvas bem distribuídas no ano. A geologia regional é constituída por uma sequência vulcânica de rochas ácidas da Formação Serra Geral, com predomínio de riodacito.
Na área de estudo foam identificados Neossolos Litólicos e Neossolos Regolíticos, Cambissolos Háplicos e Cambissolos Húmicos, Latossolos Vermelhos e Gleissolos Melânicos, os quais representam a topossequência clássica da área de estudo (Figura 1c). Conforme as tendências observadas no campo, em locais de relevo plano ou suavemente ondulado com boa drenagem estão solos profundos com sequência de horizontes A-Bw (Latossolos), em condições de má drenagem, ocorrem solos com a sequência de horizonte A-Cg (Gleissolos). Nos relevos ondulado ou fortemente ondulado, predominam solos rasos com a sequência de horizontes A ou A-Bi (Neossolos e Cambissolos).
Figura 1. Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).
Amostragem
Definimos uma malha amostral contendo 102 pontos para coleta de dados (Figura 1 b). Os pontos foram alocados pelo algoritmo Hipercubo Latino Condicionado (Conditioned Latin Hipercube Sampling – cLHS) através da função \(clhs()\) implementado no pacote clhs. Foram consideradas como variáveis condicionantes da amostragem ELEV, VD, TWI, CNBL e DECLI as quais, juntas, explicaram aproximadamente 86% da variância topográfica, identificada através da Análise de Componentes Principais (ACP).
Em cada ponto amostral foi medida a profundidade do solum (PF) e a altura das árvores (h). Consideramos solum a espessura máxima do solo onde as raízes podem se desenvolver sem impedimentos físicos para penetração livre das raízes. Os fatores limitantes considerados foram o lençol freático elevado e o contato com rocha consolidada (contato lítico) com ou sem fissuras. Para representar h utilizamos o valor médio dos quatro indivíduos de \(Pinus taeda\) mais próximos de cada ponto amostral foram mensurados. Portanto, cada ponto de amostragem representa uma área (bloco) de aproximadamente 10 x 10 metros no campo.
pontos <- read.csv('../data/GateadosDados.csv', dec = ".", sep= ";", stringsAsFactors = FALSE)
sp :: coordinates(pontos) <- c('X' , 'Y')
sp :: proj4string(pontos) <- wgs84utm22s #coordenada referencia
pontos <- sp :: spTransform(pontos, wgs84utm22s) #transforma coordenada para wgs84
pontos@data## ID CLASSE AREIA_1 ARGILA_1 h CAP PF PFd EA
## 1 1 RL 142 544 31.1 161.8 20 2.0 20
## 2 2 RL 159 432 32.8 151.3 10 1.0 20
## 3 3 RL 75 639 28.8 138.5 20 2.0 20
## 4 4 CX 121 628 29.5 162.3 100 10.0 30
## 5 5 RL 181 596 34.1 173.0 30 3.0 30
## 6 6 RL 121 647 32.9 146.8 30 3.0 30
## 7 7 RL 87 591 29.2 146.8 20 2.0 20
## 8 8 GX 93 564 35.2 159.3 100 10.0 30
## 9 9 RL 164 476 31.9 137.8 10 1.0 10
## 10 10 CX 135 558 33.1 139.3 100 10.0 50
## 11 12 RR 107 559 31.0 155.0 55 5.5 55
## 12 14 RL 178 539 29.6 136.8 10 1.0 10
## 13 15 RL 111 580 30.1 146.8 10 1.0 10
## 14 16 CX 100 578 34.5 148.3 100 10.0 40
## 15 17 RL 178 554 30.3 143.8 15 1.5 15
## 16 18 RL 141 573 24.0 132.8 10 1.0 10
## 17 19 RL 109 625 30.1 142.5 30 3.0 30
## 18 20 CX 137 375 34.0 155.0 100 10.0 35
## 19 21 CX 117 426 30.8 189.8 100 10.0 40
## 20 22 RL 128 470 31.7 151.8 20 2.0 20
## 21 23 CX 113 551 36.0 156.8 100 10.0 30
## 22 24 RL 134 534 31.5 140.8 15 1.5 15
## 23 25 CX 121 596 35.4 156.5 100 10.0 30
## 24 26 GX 99 423 32.9 166.0 100 10.0 40
## 25 27 RL 193 445 26.7 129.8 30 3.0 30
## 26 28 RR 137 550 34.7 147.3 50 5.0 50
## 27 29 RL 193 520 31.2 140.5 10 1.0 10
## 28 30 RR 156 551 31.3 139.0 40 4.0 40
## 29 31 CX 96 566 35.5 145.5 100 10.0 40
## 30 32 CX 151 534 34.0 146.3 100 10.0 30
## 31 33 CX 151 465 31.4 188.3 60 6.0 30
## 32 34 RL 152 504 25.5 148.5 10 1.0 10
## 33 35 CX 125 617 33.1 150.8 100 10.0 30
## 34 36 RR 119 627 28.5 142.8 30 3.0 30
## 35 37 CX 107 549 36.5 145.3 100 10.0 35
## 36 38 CX 157 555 34.2 139.3 100 10.0 30
## 37 39 CX 125 575 35.5 169.0 90 9.0 50
## 38 40 RL 147 498 29.0 149.0 15 1.5 15
## 39 41 RL/RR 107 580 32.6 149.6 30 3.0 30
## 40 42 CX 161 574 33.5 146.8 100 10.0 40
## 41 43 CX 127 614 34.1 142.8 100 10.0 30
## 42 44 RL 134 614 32.4 146.0 20 2.0 20
## 43 45 RR 73 330 29.1 127.8 40 4.0 40
## 44 46 CX 69 633 35.9 147.8 100 10.0 60
## 45 47 CX 81 664 33.7 136.5 100 10.0 40
## 46 48 RL 103 574 33.4 144.8 40 4.0 20
## 47 49 CX 90 643 32.1 146.0 100 10.0 40
## 48 51 RL 159 432 28.1 125.5 15 1.5 15
## 49 52 RL 181 596 30.5 123.8 35 3.5 20
## 50 53 RR 136 441 25.6 138.0 30 3.0 30
## 51 54 CX 93 641 30.0 148.5 100 10.0 30
## 52 56 RR 73 330 27.9 135.5 40 4.0 30
## 53 57 CX 40 614 36.2 171.5 100 10.0 40
## 54 58 RL 167 522 32.8 143.0 20 2.0 20
## 55 59 CX 59 432 28.8 137.8 100 10.0 40
## 56 60 RL 163 531 30.4 150.3 10 1.0 10
## 57 61 CX 51 453 32.5 176.5 100 10.0 50
## 58 62 RR 80 414 28.1 139.8 70 7.0 70
## 59 63 RL 202 512 31.5 147.0 30 3.0 20
## 60 64 CX 157 555 34.4 159.8 100 10.0 60
## 61 65 GX 127 630 30.6 152.3 70 7.0 30
## 62 66 GX 96 210 30.1 162.3 60 6.0 20
## 63 67 RR 100 471 33.1 154.0 50 5.0 50
## 64 68 RR 137 659 30.0 146.8 40 4.0 20
## 65 69 RL 229 495 24.4 140.0 10 1.0 10
## 66 70 RR 164 478 28.9 137.5 55 5.5 55
## 67 71 RR 120 633 30.0 170.3 60 6.0 20
## 68 72 CX 89 560 30.4 148.5 60 6.0 45
## 69 73 RL 152 563 27.2 160.0 10 1.0 10
## 70 74 RL 186 484 28.9 135.3 15 1.5 15
## 71 75 RL 149 549 32.2 130.0 30 3.0 30
## 72 76 RR 70 648 28.0 134.5 50 5.0 30
## 73 79 CX 163 482 30.8 151.3 60 6.0 20
## 74 80 CX 175 551 30.0 133.8 70 7.0 30
## 75 81 CX 58 499 28.9 175.0 100 10.0 40
## 76 82 CX 90 572 32.4 136.8 100 10.0 30
## 77 83 RL 75 639 28.0 130.0 20 2.0 20
## 78 84 RL 219 491 28.4 135.3 15 1.5 10
## 79 86 GX 62 425 28.8 154.5 100 10.0 40
## 80 87 CX 87 599 31.6 137.8 60 6.0 30
## 81 88 CX 117 533 31.6 132.0 70 7.0 30
## 82 89 CX 93 379 31.5 148.5 10 1.0 60
## 83 90 CX 68 600 31.5 137.3 100 10.0 40
## 84 91 RL 87 591 29.8 126.5 20 2.0 20
## 85 92 CX 92 576 28.5 160.0 100 10.0 40
## 86 93 CX 68 657 30.5 142.3 70 7.0 30
## 87 94 CX 109 611 29.8 139.5 100 10.0 40
## 88 95 CX 84 556 26.7 152.8 90 9.0 38
## 89 97 CX 84 556 31.0 140.3 80 8.0 40
## 90 99 OX 90 537 29.4 119.0 100 10.0 50
## 91 110 CX 100 558 30.8 153.3 60 6.0 20
## 92 114 RL 160 502 25.6 152.0 20 2.0 20
## 93 115 CX 80 549 29.1 134.8 90 9.0 60
## 94 122 CX 74 635 30.7 148.8 70 7.0 60
## 95 131 CX 30 714 29.2 146.8 100 10.0 40
## 96 132 CX 154 576 33.5 150.5 100 10.0 30Verificar a normalidade dos dados
stats::shapiro.test(pontos$PFd)##
## Shapiro-Wilk normality test
##
## data: pontos$PFd
## W = 0.83462, p-value = 5.627e-09stats::shapiro.test(pontos$h)##
## Shapiro-Wilk normality test
##
## data: pontos$h
## W = 0.98795, p-value = 0.5348PF bimodal h normal
par(mfrow=c(1,2))
ex <- graphics::hist(pontos$PFd, breaks=10, xlab="", col="grey90", cex.axis=1,
xlim=c(1,10), ylim=c(0,37),ylab="",main=" ")
xfit<-base::seq(min(pontos$PFd),max(pontos$PFd), length=50)
yfit<-stats::dnorm(xfit,mean=mean(pontos$PFd),sd=sd(pontos$PFd))
yfit<- yfit*diff(ex$mids[1:2])*length(pontos$PFd)
lines(xfit, yfit, col="red", lwd=2)
mtext("Frequência",line=2.6, side=2, cex=1)
mtext("PF (dm)",line=2.6, side=1, cex=1)
graphics::rug((pontos$PFd), col="red")
box()
ex2 <- graphics::hist(pontos$h, breaks=10, xlab="", col="grey90", cex.axis=1,
xlim=c(23,38), ylim=c(0,16),ylab="",main="")
xfit<-base::seq(min(pontos$h),max(pontos$h), length=50)
yfit<-stats::dnorm(xfit,mean=mean(pontos$h),sd=sd(pontos$h))
yfit<- yfit*diff(ex2$mids[1:2])*length(pontos$h)
lines(xfit, yfit, col="red", lwd=2)
mtext("Frequência",line=2.6, side=2, cex=1)
mtext("h (m)",line=2.6, side=1, cex=1)
graphics::rug((pontos$h), col="red")
box()
Utilizei o MDE disponibilizado pelo Governo do Estado de SC - Secretaria de Estado do Desenvolvimento Econômico Sustentável, proveniente do Levantamento Aerofotogramétrico em 2010. Os dados, disponibilizados com resolução de 1 metro, foram reamostrados para resolução espacial de 10 metros utilizando a ferramenta \(reamostragem\) no software SAGA GIS.
A partir do MDE foram derivadas 12 variáveis topográficas utilizando a ferramenta \(Tarrain Analyses\) do software SAGA GIS.
ASP <- raster::raster("../data/Covars/ASP.tif")
CNBL <- raster::raster("../data/Covars/CNBL.tif")
DECLI <- raster::raster("../data/Covars/DECLI.tif")
ELEV <- raster::raster("../data/Covars/ELEV.tif")
LS <- raster::raster("../data/Covars/LS.tif")
PLC <- raster::raster("../data/Covars/PLC.tif")
RSP <- raster::raster("../data/Covars/RSP.tif")
TPI <- raster::raster("../data/Covars/TPI.tif")
TRI <- raster::raster("../data/Covars/TRI.tif")
TWI <- raster::raster("../data/Covars/TWI.tif")
VD <- raster::raster("../data/Covars/VD.tif")
VDCN <- raster::raster("../data/Covars/VDCN.tif")
RFPF <- raster::raster("../data/Covars/predictionPF.tif")A partir da função extract implementada no pacote raster extraí os valores de cada objeto raster nas localizações de cada observação contidas no objeto espacial “pontos”.
pontos$ASP <- raster::extract(ASP, pontos)
pontos$CNBL <- raster::extract(CNBL, pontos)
pontos$DECLI <- raster::extract(DECLI, pontos)
pontos$ELEV <- raster::extract(ELEV, pontos)
pontos$LS <- raster::extract(LS, pontos)
pontos$PLC <- raster::extract(PLC, pontos)
pontos$RSP <- raster::extract(RSP, pontos)
pontos$TPI <- raster::extract(TPI, pontos)
pontos$TRI <- raster::extract(TRI, pontos)
pontos$TWI <- raster::extract(TWI, pontos)
pontos$VD <- raster::extract(VD, pontos)
pontos$VDCN <- raster::extract(VDCN, pontos)
pontos$RFPF <- raster::extract(RFPF, pontos)Além dos atributos do terreno, foi utilizado o mapa predito via random forest (RF) utilizando o comando \(rf()\) implementado no pacote caret, com número padrão de preditores a serem selecionados em cada nó padrão (mtry = padrão) e 1000 árvores \((ntrees = 1000)\), utilizando como preditoras as variáveis topográficas oriundas do MDE
sp::spplot(RFPF, scales = list(draw = TRUE))
Para fins de ordenamento de produção, a área possui 12 parcelas fixas de inventário contínuo (PIC) de 500 metros quadrados cada que são utilizados na estimativa da produtividade local.
Cada PIC é classificada em função da média de altura das 100 árvores de maior perímetro basal da parcela, expressa como altura dominante (Hdom). Em função da Hdom e seus incrimentos anuais é estabelecido o índices de sítio (IS). Esse índice é então atribuido à todo o talhão.
As parcelas da área foram classificação em 4 níveis, em que 1 corresponde ao sítio com melhor produtividade (Figura 2).
Figura 2 - Área de estudo, com indicação da área produtiva e parcelas de inventário contínuo contento os índices de sítio e respectivos valores de altura dominante
Os níveis de produtividade foram atribuidos aos polígonos e, naquelas em que não existe parcela de PIC, o valor de sítio foi estimado com base na altura das árvores amostradas á campo e das PCIs vizinhas.
pistola <-
raster::shapefile('../data/contorno/limite_projeto.shp',
stringsAsFactors = FALSE, encoding = 'UTF-8') %>%
sp::spTransform(wgs84utm22s)Utilizei a função \(shapefile\) implementada no pacote raster para carregar o polígono da área de estudo - armazenado no objeto pistola e o polígono com informações das sítio - armazenado no objeto ProdutividadePistola. A função sp::spTransform foi usada para projetar as coordenadas original no plano cartesiano (UTM).
ProdutividadePistola <-
raster::shapefile('../data/Produtividade/ProdutividadePistola.shp') %>%
sp::spTransform(wgs84utm22s)
sp::spplot(ProdutividadePistola, main = "Mapa de índices de Sítio")
#default.stringsAsFactors()
#str(ProdutividadePistola)Usei a função \(over\) implementada no pacote sp para identificar o nível de produtividade do polígono dentro da qual cada observação se encontra e O resultado foi armazenado em uma coluna definida como Sitio no objeto mesmo objeto espacial pontos.
pontos$Sitio <- sp::over(x = pontos, y = ProdutividadePistola, na.omit = TRUE) %>% unlist()
pontos$Sitio <- as.factor(pontos@data$Sitio)col_sitio <- terrain.colors(nlevels(ProdutividadePistola$ProdutividadePistola))
sp::spplot(
ProdutividadePistola, scales = list(draw = T),
main = 'Localização dos pontos') +
lattice::xyplot(Y ~ X, data = as.data.frame(pontos@coords),
pch = 20, col = 'red', lwd = 2, cex = 1.5) %>%
latticeExtra::as.layer()
Parte I - Profundidade do solo
n = 96 observações para calibração + predição + validação (108 ha) Validação: validação cruzada (leave-one-out)
Modelo linear misto de variação espacial (georob)
Para utilizar esse modelo, supus que os dados são uma realização de um campo aleatório com distribuição normal que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos aleatórios e erro aleatório independente. Assim, considerei como efeitos fixos da variação da PFd a predição da PF via floresta aleatória, atributos de terreno oriundos do MDE e os índices de sítio.
Ajuste do variograma:
limites <- seq(0, 2000, length.out = 15)
vario<- georob::sample.variogram(PFd ~ RFPF + ELEV + DECLI + Sitio,
data= pontos, locations = ~ X + Y, lag.dist.def = limites, estimator = c( "matheron"),
mean.angle = TRUE, annotate.npairs = T) %>%
plot()
lags <- seq(0, 2000, length.out = 15)
georob::sample.variogram(
PFd ~ RFPF + ELEV + DECLI, data = pontos, locations = ~ X + Y, lag.dist.def = lags,
xy.angle.def = c(0, 22.5, 67.5, 112.5, 157.5, 180)) %>%
plot(type = "b", ylab = 'Semivariância', xlab = 'Distância de separação (km)')
Avaliada a evolução da semivariância nas diferentes direções, há evidência da existência de estruturas de autocorrelação espacial dependentes da direção, então o processo espacial é anisotrópico. Porém, assumi a isotropia do processo espacial e o semivariograma amostral gerado foi independente da direção.
Ajustei ao variograma amostral um modelo exponencial do variograma usando o método dos quadrados mínimos não-lineares ponderados, com ponderação definida conforme o método de Cressie. O processo de estimativa dos parâmetros do modelo exponencial do variograma foi conduzido via otimização usando a função stats::optim(method = “BFGS”).
produzidos pelo otimizador a cada 10 iterações:
vario_fit <-
georob::fit.variogram.model(
vario, variogram.model = 'RMexp', param = c(variance = 13, nugget = 10, scale = 70),
weighting.method = "cressie", method = "BFGS")##
##
## variance snugget nugget scale
## Variogram parameters 0.000000e+00 0.000000e+00 0.000000e+00 3.204568e+22
##
##
## variance snugget nugget scale
## Variogram parameters 2.731246e+28 0.000000e+00 2.687536e-41 5.213753e+02summary(vario_fit)##
## Call:georob::fit.variogram.model(sv = vario, variogram.model = "RMexp",
## param = c(variance = 13, nugget = 10, scale = 70), weighting.method = "cressie",
## method = "BFGS")
##
## Convergence in 212 function and 80 Jacobian/gradient evaluations
##
## Residual Sum of Squares: 41.57611
##
## Residuals (epsilon):
## Min 1Q Median 3Q Max
## -2.37531 -1.66544 -0.61869 0.08033 1.99610
##
## Variogram: RMexp
## Estimate Lower Upper
## variance 1.312e+01 1.283e+01 1.341e+01
## snugget(fixed) 0.000e+00 NA NA
## nugget 5.055e-04 1.088e-46 2.349e+39
## scale 6.635e+01 5.482e+01 8.029e+01plot(vario, type = "b", xlab = 'Distância de separação (m)', ylab = 'Semivariância')
lines(vario_fit, col = "red", lty = 'dashed')
Na geoestatística clássica, a variância não explicada é expressa em um único parâmetro nugget. Porém, no modelo linear misto podemos separar o parâmetro nugget em dois componentes: a variãncia devida aos erros de medida, modelada pelo parâmetro nugget e variância devida à variação espacial em pequena escala, modelada pelo parâmetro snugget. Considerando que durante a obteção dos dados de campo as informações de PF os valores foram arredondadas pelas em decímetros, considerei que a variância do erro de medida é igual a 0.5, ou seja, 62% do parâmetro nugget. A variância restante foi atribuída a variação espacial não auto-correlacionada espacialmente - não capturada pelo plano amostral snugget.
Feitas essa suposição, ajustei novamente o modelo novamente, agora mantendo ambos nugget e snugget fixos.
vario_fit_error <- georob::georob(
PFd ~ RFPF, pontos, locations = ~ X + Y, variogram.model = 'RMexp',
param = c(variance = vario_fit$variogram.object[[1]]$param[['variance']],
nugget = vario_fit$variogram.object[[1]]$param[['nugget']] * 0.5,
snugget = vario_fit$variogram.object[[1]]$param[['nugget']] * 0.5,
scale = vario_fit$variogram.object[[1]]$param[['scale']]),
fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE),
tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(vario_fit_error)##
## Call:georob::georob(formula = PFd ~ RFPF, data = pontos, locations = ~X +
## Y, variogram.model = "RMexp", param = c(variance = vario_fit$variogram.object[[1]]$param[["variance"]],
## nugget = vario_fit$variogram.object[[1]]$param[["nugget"]] *
## 0.5, snugget = vario_fit$variogram.object[[1]]$param[["nugget"]] *
## 0.5, scale = vario_fit$variogram.object[[1]]$param[["scale"]]),
## fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE),
## tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
##
## Tuning constant: 1000
##
## Convergence in 7 function and 6 Jacobian/gradient evaluations
##
## Estimating equations (gradient)
##
## xi scale
## Gradient : -1.073960e-02 1.324515e-02
##
## Maximized restricted log-likelihood: -205.2834
##
## Predicted latent variable (B):
## Min 1Q Median 3Q Max
## -6.86420 -0.67666 0.09004 1.06093 6.73776
##
## Residuals (epsilon):
## Min 1Q Median 3Q Max
## -4.251e-04 -5.509e-05 7.842e-06 6.590e-05 4.148e-04
##
## Standardized residuals:
## Min 1Q Median 3Q Max
## -3.4088 -0.4274 0.0635 0.5154 3.3177
##
##
## Gaussian REML estimates
##
## Variogram: RMexp
## Estimate Lower Upper
## variance 4.172e+00 3.103e+00 5.608
## snugget(fixed) 2.527e-04 NA NA
## nugget(fixed) 2.527e-04 NA NA
## scale 3.390e+01 1.829e+01 62.863
##
##
## Fixed effects coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.64442 0.73305 -4.972 2.98e-06 ***
## RFPF 0.15278 0.01117 13.674 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error (sqrt(nugget)): 0.0159
##
## Robustness weights:
## 46 weights are ~= 1. The remaining 50 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.9423 0.9889 0.9956 0.9902 0.9978 0.9990plot(vario, type = "b", xlab = 'Distância de separação (m)', ylab = 'Semivariância', lty = 'dashed')
lines(vario_fit, col = "navyblue")
lines(vario_fit_error, col = "orange")
CRIAR O GRID: Suporte de predição pontual
grid <- sp::spsample(ProdutividadePistola, 10000, type = 'regular')
colnames(grid@coords) <- colnames(pontos@coords)grid$Sitio <- sp::over(grid, ProdutividadePistola)
grid$DECLI <- raster::extract(DECLI, grid)
grid$ELEV <- raster::extract(ELEV, grid)
grid$RFPF <- raster::extract(RFPF, grid)
grid## class : SpatialPointsDataFrame
## features : 10004
## extent : 516813.9, 518852.5, 6908631, 6909820 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0
## variables : 4
## names : Sitio, DECLI, ELEV, RFPF
## min values : 1, 0.000142422635690309, 873.203247070312, 16.9653333333333
## max values : 4, 0.459370046854019, 936.466369628906, 94.5196666666667grid <-
sp::SpatialPointsDataFrame(
coords = grid@coords,
data = data.frame(grid),
proj4string = grid@proj4string)
colnames(grid@coords) <- colnames(pontos@coords)Argumento type: fazer a predição do sinal “signal” z, “response” se for y ou efeitos fixos “trend” se quiser predizer erro de medida, y = z + e, então response - signal = erro se tendência - predito = componente aleatório
Usar control.predict.georob(extended.output = TRUE)
pred_ponto <- predict(
vario_fit_error, newdata = grid, type= "response", signif = 0.95,
control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE## Warning in points2grid(points, tolerance, round): grid has empty column/
## rows in dimension 1str(pred_ponto)## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
## ..@ data :'data.frame': 10004 obs. of 8 variables:
## .. ..$ pred : num [1:10004] 7.87 8.25 3.92 4.24 8.3 ...
## .. ..$ se : num [1:10004] 2.06 2.06 2.03 2.02 2.06 ...
## .. ..$ lower : num [1:10004] 3.8387 4.2183 -0.0577 0.2717 4.2687 ...
## .. ..$ upper : num [1:10004] 11.9 12.3 7.9 8.2 12.3 ...
## .. ..$ trend : num [1:10004] 7.94 8.31 4.08 4.4 8.4 ...
## .. ..$ var.pred : num [1:10004] 0.101 0.113 0.171 0.186 0.122 ...
## .. ..$ cov.pred.target: num [1:10004] 0.0231 0.0263 0.1114 0.1332 0.0347 ...
## .. ..$ var.target : num [1:10004] 4.17 4.17 4.17 4.17 4.17 ...
## .. ..- attr(*, "variogram.object")=List of 1
## .. .. ..$ :List of 9
## .. .. .. ..$ variogram.model: chr "RMexp"
## .. .. .. ..$ param : Named num [1:4] 4.17 2.53e-04 2.53e-04 3.39e+01
## .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
## .. .. .. ..$ fit.param : Named logi [1:4] TRUE FALSE FALSE TRUE
## .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
## .. .. .. ..$ isotropic : logi TRUE
## .. .. .. ..$ aniso : Named num [1:5] 1 1 90 90 0
## .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
## .. .. .. ..$ fit.aniso : Named logi [1:5] FALSE FALSE FALSE FALSE FALSE
## .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
## .. .. .. ..$ sincos :List of 6
## .. .. .. .. ..$ co: num 6.12e-17
## .. .. .. .. ..$ so: num 1
## .. .. .. .. ..$ cp: num 6.12e-17
## .. .. .. .. ..$ sp: num 1
## .. .. .. .. ..$ cz: num 1
## .. .. .. .. ..$ sz: num 0
## .. .. .. ..$ rotmat : num [1:2, 1:2] 1.00 -6.12e-17 6.12e-17 1.00
## .. .. .. ..$ sclmat : Named num [1:2] 1 1
## .. .. .. .. ..- attr(*, "names")= chr [1:2] "" "f1"
## .. ..- attr(*, "psi.func")= chr "logistic"
## .. ..- attr(*, "tuning.psi")= num 1000
## .. ..- attr(*, "type")= chr "response"
## ..@ coords.nrs : num(0)
## ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots
## .. .. ..@ cellcentre.offset: Named num [1:2] 516814 6908631
## .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
## .. .. ..@ cellsize : Named num [1:2] 8.68 8.68
## .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
## .. .. ..@ cells.dim : Named int [1:2] 236 138
## .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
## ..@ grid.index : int [1:10004] 32416 32417 32139 32140 32180 32181 31903 31904 31905 31943 ...
## ..@ coords : num [1:10004, 1:2] 517534 517543 517178 517187 517534 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:10004] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "X" "Y"
## ..@ bbox : num [1:2, 1:2] 516810 6908627 518857 6909824
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "X" "Y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NAspplot(pred_ponto)
at <- pred_ponto@data[, c("pred", "lower", "upper")] %>% range()
sp::spplot(pred_ponto, zcol = c("lower", "pred", "upper"), main = "prediction")
sp::spplot(pred_ponto, zcol = 'se', main = "SE")
validacao <- georob::cv(vario_fit_error, nset = 95)
summary(validacao)##
## Statistics of cross-validation prediction errors
## me mede rmse made qne msse medsse
## 0.001978 0.129553 2.008528 1.379093 1.583218 1.072934 0.229126
## crps
## 1.083096